#### Setting environment ####
require(estimatr)
require(dichromat)
color <- paste0(colorschemes$Categorical.12, "80")

# Function to compute summary statistics for a given variable
summary.statistics <- function(x) {
  out <- rep(NA, 5)
  names(out) <- c("Mean", "SD", "Min", "Max", "N")
  out[1] <- mean(x, na.rm = TRUE)
  out[2] <- sd(x, na.rm = TRUE)
  out[3] <- min(x, na.rm = TRUE)
  out[4] <- max(x, na.rm = TRUE)
  out[5] <- sum(! is.na(x))
  round(out, 2)
}

# Function to compute conditional average treatment effects (CATEs)
conditional.effects <- function(results, treatment, covariate, cov.values) {
  int.term <- paste0(treatment, ":", covariate)
  est <- results$coefficients[treatment] + 
    results$coefficients[int.term] * cov.values
  vcv <- vcov(results)
  se <- sqrt(vcv[treatment, treatment] + 
               (cov.values ^ 2) * vcv[int.term, int.term] + 
               2 * cov.values * vcv[treatment, int.term])
  lwr <- est + qnorm(0.025) * se
  upr <- est + qnorm(0.975) * se
  cbind(est, lwr, upr)
}

# Function to compute CATEs while holding another variable at a fixed value
conditional.effects.fixed <- function(results, treatment, covariate, cov.values, 
                                      fixed.var, fixed.val) {
  int.term <- paste0(treatment, ":", covariate)
  int.term.fixed <- paste0(treatment, ":", fixed.var)
  est <- results$coefficients[treatment] + 
    results$coefficients[int.term] * cov.values + 
    results$coefficients[int.term.fixed] * fixed.val
  vcv <- vcov(results)
  se <- sqrt(vcv[treatment, treatment] + 
               (cov.values ^ 2) * vcv[int.term, int.term] + 
               (fixed.val ^ 2) * vcv[int.term.fixed, int.term.fixed] + 
               2 * cov.values * vcv[treatment, int.term] + 
               2 * fixed.val * vcv[treatment, int.term.fixed] + 
               2 * cov.values * fixed.val * vcv[int.term, int.term.fixed])
  lwr <- est + qnorm(0.025) * se
  upr <- est + qnorm(0.975) * se
  cbind(est, lwr, upr)
}

#### Data ####
data <- read.csv("replication_data.csv")

# Number of respondents
nrow(data)

#### Descriptive statistics ####
## Summary statistics (Table A1)
# Respondent attributes
summary.statistics(data$R.woman)
summary.statistics(data$R.age)
summary.statistics(data$tenure)
summary.statistics(data$DID)
summary.statistics(data$magnitude)

# Outcome variables
summary.statistics(data$nomination)
summary.statistics(data$support)
summary.statistics(data$electability)

## Representativeness comparison (Figure A1)
hist.gender.sample <- prop.table(table(data$R.woman, useNA = "always"))
# Population data source:
# https://www.si-gichokai.jp/research/zokusei/__icsFiles/afieldfile/2022/11/28/R04zokuseisirabe_1.pdf
hist.gender.sample.pop <- c(15293, 3307, 0) / 18600

# First panel of Figure A1
# (Remaining panels are not reproducible due to anonymity restrictions)
png("Figure_A1.png", width = 2, height = 2, units = "in", pointsize = 8, res = 2000)
par(mar = c(2, 3.5, 2, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, length(hist.gender.sample) + 0.5), 
     ylim = c(0, max(c(hist.gender.sample, hist.gender.sample.pop))), 
     main = "Gender", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:length(hist.gender.sample)) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(0, hist.gender.sample[i], hist.gender.sample[i], 0), 
          border = color[12], col = color[11])
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(0, hist.gender.sample.pop[i], hist.gender.sample.pop[i], 0), 
          density = 30, border = color[10], col = color[9])
}
mtext(c("Man", "Woman", "NA"), 1, at = 1:length(hist.gender.sample), cex = 0.7)
mtext("Density", 2, 2.5, cex = 0.7)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### H1 and H2 results ####
# Table A2
nomination.base.results <- lm_robust(nomination ~ woman + young, 
                                     data, clusters = code, 
                                     fixed_effects = cluster)
nomination.base.summary <- summary(nomination.base.results)
round(nomination.base.summary$coefficients[, c(1, 2, 4)], 3)
nomination.base.results$nobs

support.base.results <- lm_robust(support ~ woman + young, 
                                  data, clusters = code, 
                                  fixed_effects = cluster)
support.base.summary <- summary(support.base.results)
round(support.base.summary$coefficients[, c(1, 2, 4)], 3)
support.base.results$nobs

electability.base.results <- lm_robust(electability ~ woman + young, 
                                       data, clusters = code, 
                                       fixed_effects = cluster)
electability.base.summary <- summary(electability.base.results)
round(electability.base.summary$coefficients[, c(1, 2, 4)], 3)
electability.base.results$nobs

# Figure 1
png("Figure_1.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(0, 0.25), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(0.05, 0.25, 0.05), lty = 3, col = "lightgray")
segments(1, nomination.base.results$conf.low[1], 
         1, nomination.base.results$conf.high[1])
points(1, nomination.base.results$coefficients[1], pch = 19)
segments(2, support.base.results$conf.low[1], 
         2, support.base.results$conf.high[1])
points(2, support.base.results$coefficients[1], pch = 19)
segments(3, electability.base.results$conf.low[1], 
         3, electability.base.results$conf.high[1])
points(3, electability.base.results$coefficients[1], pch = 19)
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.base.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.base.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.base.summary$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### H3 results ####
# Table A3
nomination.int.results <- lm_robust(nomination ~ woman * young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
nomination.int.summary <- summary(nomination.int.results)
round(nomination.int.summary$coefficients[, c(1, 2, 4)], 3)
nomination.int.results$nobs

support.int.results <- lm_robust(support ~ woman * young, 
                                 data, clusters = code, 
                                 fixed_effects = cluster)
support.int.summary <- summary(support.int.results)
round(support.int.summary$coefficients[, c(1, 2, 4)], 3)
support.int.results$nobs

electability.int.results <- lm_robust(electability ~ woman * young, 
                                      data, clusters = code, 
                                      fixed_effects = cluster)
electability.int.summary <- summary(electability.int.results)
round(electability.int.summary$coefficients[, c(1, 2, 4)], 3)
electability.int.results$nobs

# Regressions with the young dummy reversed (for Figure 2)
nomination.int.results.rev <- lm_robust(nomination ~ woman * I(1 - young), 
                                        data, clusters = code, 
                                        fixed_effects = cluster)

support.int.results.rev <- lm_robust(support ~ woman * I(1 - young), 
                                     data, clusters = code, 
                                     fixed_effects = cluster)

electability.int.results.rev <- lm_robust(electability ~ woman * I(1 - young), 
                                          data, clusters = code, 
                                          fixed_effects = cluster)

# Figure 2
png("Figure_2.png",width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(-0.1, 0.3), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.1, 0.3, 0.1)[-2], lty = 3, col = "lightgray")
segments(0.85, nomination.int.results.rev$conf.low[1], 
         0.85, nomination.int.results.rev$conf.high[1])
points(0.85, nomination.int.results.rev$coefficients[1], pch = 19)
text(0.85, nomination.int.results.rev$conf.high[1], "Younger", pos = 3)
segments(1.15, nomination.int.results$conf.low[1], 
         1.15, nomination.int.results$conf.high[1])
points(1.15, nomination.int.results$coefficients[1], pch = 21, bg = "white")
text(1.15, nomination.int.results$conf.low[1], "Older", pos = 1)
segments(1.85, support.int.results.rev$conf.low[1], 
         1.85, support.int.results.rev$conf.high[1])
points(1.85, support.int.results.rev$coefficients[1], pch = 19)
segments(2.15, support.int.results$conf.low[1], 
         2.15, support.int.results$conf.high[1])
points(2.15, support.int.results$coefficients[1], pch = 21, bg = "white")
segments(2.85, electability.int.results.rev$conf.low[1], 
         2.85, electability.int.results.rev$conf.high[1])
points(2.85, electability.int.results.rev$coefficients[1], pch = 19)
segments(3.15, electability.int.results$conf.low[1], 
         3.15, electability.int.results$conf.high[1])
points(3.15, electability.int.results$coefficients[1], pch = 21, bg = "white")
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.int.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.int.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.int.summary$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Heterogeneity analysis by respondent gender ####
# Table A4
nomination.gender.results <- lm_robust(nomination ~ woman * R.woman + young, 
                                       data, clusters = code, 
                                       fixed_effects = cluster)
nomination.gender.summary <- summary(nomination.gender.results)
round(nomination.gender.summary$coefficients[, c(1, 2, 4)], 3)
nomination.gender.results$nobs

support.gender.results <- lm_robust(support ~ woman * R.woman + young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
support.gender.summary <- summary(support.gender.results)
round(support.gender.summary$coefficients[, c(1, 2, 4)], 3)
support.gender.results$nobs

electability.gender.results <- lm_robust(electability ~ woman * R.woman + young, 
                                         data, clusters = code, 
                                         fixed_effects = cluster)
electability.gender.summary <- summary(electability.gender.results)
round(electability.gender.summary$coefficients[, c(1, 2, 4)], 3)
electability.gender.results$nobs

# Regressions with respondent gender dummy reversed (for Figure 3)
nomination.gender.results.rev <- lm_robust(nomination ~ woman * I(1 - R.woman) + young, 
                                           data, clusters = code, 
                                           fixed_effects = cluster)

support.gender.results.rev <- lm_robust(support ~ woman * I(1 - R.woman) + young, 
                                        data, clusters = code, 
                                        fixed_effects = cluster)

electability.gender.results.rev <- lm_robust(electability ~ woman * I(1 - R.woman) + young, 
                                             data, clusters = code, 
                                             fixed_effects = cluster)

# Figure 3
png("Figure_3.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(-0.1, 0.4), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.1, 0.4, 0.1)[-2], lty = 3, col = "lightgray")
segments(0.85, nomination.gender.results.rev$conf.low[1], 
         0.85, nomination.gender.results.rev$conf.high[1])
points(0.85, nomination.gender.results.rev$coefficients[1], pch = 19)
text(0.85, nomination.gender.results.rev$conf.high[1], "Female", pos = 3)
segments(1.15, nomination.gender.results$conf.low[1], 
         1.15, nomination.gender.results$conf.high[1])
points(1.15, nomination.gender.results$coefficients[1], pch = 21, bg = "white")
text(1.15, nomination.gender.results$conf.low[1], "Male", pos = 1)
segments(1.85, support.gender.results.rev$conf.low[1], 
         1.85, support.gender.results.rev$conf.high[1])
points(1.85, support.gender.results.rev$coefficients[1], pch = 19)
segments(2.15, support.gender.results$conf.low[1], 
         2.15, support.gender.results$conf.high[1])
points(2.15, support.gender.results$coefficients[1], pch = 21, bg = "white")
segments(2.85, electability.gender.results.rev$conf.low[1], 
         2.85, electability.gender.results.rev$conf.high[1])
points(2.85, electability.gender.results.rev$coefficients[1], pch = 19)
segments(3.15, electability.gender.results$conf.low[1], 
         3.15, electability.gender.results$conf.high[1])
points(3.15, electability.gender.results$coefficients[1], pch = 21, bg = "white")
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.gender.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.gender.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.gender.summary$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Heterogeneity analysis by respondent age ####
# Note: These analyses use coarsened variables and do not replicate Table A5.
nomination.age.results <- lm_robust(nomination ~ woman * R.age + young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
nomination.age.summary <- summary(nomination.age.results)
round(nomination.age.summary$coefficients[, c(1, 2, 4)], 3)

support.age.results <- lm_robust(support ~ woman * R.age + young, 
                                 data, clusters = code, 
                                 fixed_effects = cluster)
support.age.summary <- summary(support.age.results)
round(support.age.summary$coefficients[, c(1, 2, 4)], 3)

electability.age.results <- lm_robust(electability ~ woman * R.age + young, 
                                      data, clusters = code, 
                                      fixed_effects = cluster)
electability.age.summary <- summary(electability.age.results)
round(electability.age.summary$coefficients[, c(1, 2, 4)], 3)

# Compute conditional ATEs by age group
nomination.age.plot <- conditional.effects(nomination.age.results, 
                                           "woman", "R.age", 1:5)
support.age.plot <- conditional.effects(support.age.results, 
                                        "woman", "R.age", 1:5)
electability.age.plot <- conditional.effects(electability.age.results, 
                                             "woman", "R.age", 1:5)

# Create background histogram for Figure A9
hist.age <- prop.table(table(data$R.age))

# Figure A9 (uses coarsened variable; not a reproduction of Figure 4)
png("Figure_A9.png", width = 6, height = 2, units = "in", pointsize = 8, res = 2000)
layout(matrix(1:3, 1, 3))
par(mar = c(3.5, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 5.5), ylim = c(-0.2, 0.6), 
     main = "Nomination willingness", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.6, 0.2)[-2], lty = 3, col = "lightgray")
for (i in 1:5) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(-0.2, hist.age[i] / 2 - 0.2, 
            hist.age[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:5, nomination.age.plot[, 2], 1:5, nomination.age.plot[, 3])
points(1:5, nomination.age.plot[, 1], pch = 19)
text(5, 0.56, 
     paste0("N = ", format(nomination.age.summary$nobs, big.mark = ",")))
mtext("Respondent age", 1, 2.5, cex = 0.8)
axis(1, 1:5, c("20\u201330s", "40s", "50s", "60s", "70\u2013"), 
     lwd = 0.5, cex = 0.8)
mtext("Effect of female canageate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 5.5), ylim = c(-0.2, 0.6), 
     main = "Perceived support", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.6, 0.2)[-2], lty = 3, col = "lightgray")
for (i in 1:5) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(-0.2, hist.age[i] / 2 - 0.2, 
            hist.age[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:5, support.age.plot[, 2], 1:5, support.age.plot[, 3])
points(1:5, support.age.plot[, 1], pch = 19)
text(5, 0.56, 
     paste0("N = ", format(support.age.summary$nobs, big.mark = ",")))
mtext("Respondent age", 1, 2.5, cex = 0.8)
axis(1, 1:5, c("20\u201330s", "40s", "50s", "60s", "70\u2013"), 
     lwd = 0.5, cex = 0.8)
mtext("Effect of female canageate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 5.5), ylim = c(-0.2, 0.6), 
     main = "Perceived electability", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.6, 0.2)[-2], lty = 3, col = "lightgray")
for (i in 1:5) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(-0.2, hist.age[i] / 2 - 0.2, 
            hist.age[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:5, electability.age.plot[, 2], 1:5, electability.age.plot[, 3])
points(1:5, electability.age.plot[, 1], pch = 19)
text(5, 0.56, 
     paste0("N = ", format(electability.age.summary$nobs, big.mark = ",")))
mtext("Respondent age", 1, 2.5, cex = 0.8)
axis(1, 1:5, c("20\u201330s", "40s", "50s", "60s", "70\u2013"), 
     lwd = 0.5, cex = 0.8)
mtext("Effect of female canageate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Heterogeneity analysis by the urbanization score ####
# Note: These analyses use coarsened variables and do not replicate Table A6.
nomination.DID.results <- lm_robust(nomination ~ woman * (R.woman + DID) + young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
nomination.DID.summary <- summary(nomination.DID.results)
round(nomination.DID.summary$coefficients[, c(1, 2, 4)], 3)

support.DID.results <- lm_robust(support ~ woman * (R.woman + DID) + young, 
                                 data, clusters = code, 
                                 fixed_effects = cluster)
support.DID.summary <- summary(support.DID.results)
round(support.DID.summary$coefficients[, c(1, 2, 4)], 3)

electability.DID.results <- lm_robust(electability ~ woman * (R.woman + DID) + young, 
                                      data, clusters = code, 
                                      fixed_effects = cluster)
electability.DID.summary <- summary(electability.DID.results)
round(electability.DID.summary$coefficients[, c(1, 2, 4)], 3)

# Compute conditional ATEs by the urbanization score
nomination.DID.plot <- conditional.effects.fixed(nomination.DID.results, 
                                                 "woman", "DID", 1:10, 
                                                 "R.woman", 0)
support.DID.plot <- conditional.effects.fixed(support.DID.results, 
                                              "woman", "DID", 1:10, 
                                              "R.woman", 0)
electability.DID.plot <- conditional.effects.fixed(electability.DID.results, 
                                                   "woman", "DID", 1:10, 
                                                   "R.woman", 0)

# Create background histogram for Figure A10
hist.DID <- hist(data$DID, breaks = 0:10, plot = FALSE)

# Figure A10 (uses coarsened variable; not a reproduction of Figure 5)
png("Figure_A10.png", width = 6, height = 2, units = "in", pointsize = 8, res = 2000)
layout(matrix(1:3, 1, 3))
par(mar = c(3.5, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 10.5), ylim = c(-0.2, 0.4), 
     main = "Nomination willingness", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.DID$breaks) - 1)) {
  polygon(c(hist.DID$breaks[i] + 0.5, hist.DID$breaks[i] + 0.5, 
            hist.DID$breaks[i + 1] + 0.5, hist.DID$breaks[i + 1] + 0.5), 
          c(-0.2, hist.DID$density[i] - 0.2, 
            hist.DID$density[i] - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:10, nomination.DID.plot[, 2], 1:10, nomination.DID.plot[, 3]) 
points(1:10, nomination.DID.plot[, 1], pch = 19)
text(1.5, 0.37, 
     paste0("N = ", format(nomination.DID.summary$nobs, big.mark = ",")))
mtext("Urbanization score", 1, 2.5, cex = 0.8)
axis(1, 1:10, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 10.5), ylim = c(-0.2, 0.4), 
     main = "Perceived support", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.DID$breaks) - 1)) {
  polygon(c(hist.DID$breaks[i] + 0.5, hist.DID$breaks[i] + 0.5, 
            hist.DID$breaks[i + 1] + 0.5, hist.DID$breaks[i + 1] + 0.5), 
          c(-0.2, hist.DID$density[i] - 0.2, 
            hist.DID$density[i] - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:10, support.DID.plot[, 2], 1:10, support.DID.plot[, 3]) 
points(1:10, support.DID.plot[, 1], pch = 19)
text(1.5, 0.37, 
     paste0("N = ", format(support.DID.summary$nobs, big.mark = ",")))
mtext("Urbanization score", 1, 2.5, cex = 0.8)
axis(1, 1:10, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 10.5), ylim = c(-0.2, 0.4), 
     main = "Perceived electability", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.DID$breaks) - 1)) {
  polygon(c(hist.DID$breaks[i] + 0.5, hist.DID$breaks[i] + 0.5, 
            hist.DID$breaks[i + 1] + 0.5, hist.DID$breaks[i + 1] + 0.5), 
          c(-0.2, hist.DID$density[i] - 0.2, 
            hist.DID$density[i] - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:10, electability.DID.plot[, 2], 1:10, electability.DID.plot[, 3]) 
points(1:10, electability.DID.plot[, 1], pch = 19)
text(1.5, 0.37, 
     paste0("N = ", format(electability.DID.summary$nobs, big.mark = ",")))
mtext("Urbanization score", 1, 2.5, cex = 0.8)
axis(1, 1:10, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Another approach regarding the clustered robust standard errors ####
## H1 and H2
nomination.base.results.RC <- lm_robust(nomination ~ woman + young, 
                                        data, subset = code != 999, 
                                        clusters = code, fixed_effects = cluster)
nomination.base.summary.RC <- summary(nomination.base.results.RC)
round(nomination.base.summary.RC$coefficients[, c(1, 2, 4)], 3)

support.base.results.RC <- lm_robust(support ~ woman + young, 
                                     data, subset = code != 999, 
                                     clusters = code, fixed_effects = cluster)
support.base.summary.RC <- summary(support.base.results.RC)
round(support.base.summary.RC$coefficients[, c(1, 2, 4)], 3)

electability.base.results.RC <- lm_robust(electability ~ woman + young, 
                                          data, subset = code != 999,
                                          clusters = code, fixed_effects = cluster)
electability.base.summary.RC <- summary(electability.base.results.RC)
round(electability.base.summary.RC$coefficients[, c(1, 2, 4)], 3)

# Figure A2
png("Figure_A2.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(0, 0.25), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(0.05, 0.25, 0.05), lty = 3, col = "lightgray")
segments(1, nomination.base.results.RC$conf.low[1], 
         1, nomination.base.results.RC$conf.high[1])
points(1, nomination.base.results.RC$coefficients[1], pch = 19)
segments(2, support.base.results.RC$conf.low[1], 
         2, support.base.results.RC$conf.high[1])
points(2, support.base.results.RC$coefficients[1], pch = 19)
segments(3, electability.base.results.RC$conf.low[1], 
         3, electability.base.results.RC$conf.high[1])
points(3, electability.base.results.RC$coefficients[1], pch = 19)
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.base.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.base.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.base.summary.RC$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

## H3
nomination.int.results.RC <- lm_robust(nomination ~ woman * young, 
                                       data, subset = code != 999, 
                                       clusters = code, fixed_effects = cluster)
nomination.int.summary.RC <- summary(nomination.int.results.RC)
round(nomination.int.summary.RC$coefficients[, c(1, 2, 4)], 3)

support.int.results.RC <- lm_robust(support ~ woman * young, 
                                    data, subset = code != 999, 
                                    clusters = code, fixed_effects = cluster)
support.int.summary.RC <- summary(support.int.results.RC)
round(support.int.summary.RC$coefficients[, c(1, 2, 4)], 3)

electability.int.results.RC <- lm_robust(electability ~ woman * young, 
                                         data, subset = code != 999, 
                                         clusters = code, fixed_effects = cluster)
electability.int.summary.RC <- summary(electability.int.results.RC)
round(electability.int.summary.RC$coefficients[, c(1, 2, 4)], 3)

# Regressions with the young dummy reversed (for Figure A3)
nomination.int.results.RC.rev <- lm_robust(nomination ~ woman * I(1 - young), 
                                           data, subset = code != 999, 
                                           clusters = code, fixed_effects = cluster)

support.int.results.RC.rev <- lm_robust(support ~ woman * I(1 - young), 
                                        data, subset = code != 999, 
                                        clusters = code, fixed_effects = cluster)

electability.int.results.RC.rev <- lm_robust(electability ~ woman * I(1 - young), 
                                             data, subset = code != 999, 
                                             clusters = code, fixed_effects = cluster)

# Figure A3
png("Figure_A3.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(-0.1, 0.4), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.1, 0.4, 0.1)[-2], lty = 3, col = "lightgray")
segments(0.85, nomination.int.results.RC.rev$conf.low[1], 
         0.85, nomination.int.results.RC.rev$conf.high[1])
points(0.85, nomination.int.results.RC.rev$coefficients[1], pch = 19)
text(0.85, nomination.int.results.RC.rev$conf.high[1], "Younger", pos = 3)
segments(1.15, nomination.int.results.RC$conf.low[1], 
         1.15, nomination.int.results.RC$conf.high[1])
points(1.15, nomination.int.results.RC$coefficients[1], pch = 21, bg = "white")
text(1.15, nomination.int.results.RC$conf.low[1], "Older", pos = 1)
segments(1.85, support.int.results.RC.rev$conf.low[1], 
         1.85, support.int.results.RC.rev$conf.high[1])
points(1.85, support.int.results.RC.rev$coefficients[1], pch = 19)
segments(2.15, support.int.results.RC$conf.low[1], 
         2.15, support.int.results.RC$conf.high[1])
points(2.15, support.int.results.RC$coefficients[1], pch = 21, bg = "white")
segments(2.85, electability.int.results.RC.rev$conf.low[1], 
         2.85, electability.int.results.RC.rev$conf.high[1])
points(2.85, electability.int.results.RC.rev$coefficients[1], pch = 19)
segments(3.15, electability.int.results.RC$conf.low[1], 
         3.15, electability.int.results.RC$conf.high[1])
points(3.15, electability.int.results.RC$coefficients[1], pch = 21, bg = "white")
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.int.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.int.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.int.summary.RC$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

## Heterogeneity analysis by respondents' gender
nomination.gender.results.RC <- lm_robust(nomination ~ woman * R.woman + young, 
                                          data, subset = code != 999, 
                                          clusters = code, fixed_effects = cluster)
nomination.gender.summary.RC <- summary(nomination.gender.results.RC)
round(nomination.gender.summary.RC$coefficients[, c(1, 2, 4)], 3)

support.gender.results.RC <- lm_robust(support ~ woman * R.woman + young, 
                                       data, subset = code != 999, 
                                       clusters = code, fixed_effects = cluster)
support.gender.summary.RC <- summary(support.gender.results.RC)
round(support.gender.summary.RC$coefficients[, c(1, 2, 4)], 3)

electability.gender.results.RC <- lm_robust(electability ~ woman * R.woman + young, 
                                            data, subset = code != 999, 
                                            clusters = code, fixed_effects = cluster)
electability.gender.summary.RC <- summary(electability.gender.results.RC)
round(electability.gender.summary.RC$coefficients[, c(1, 2, 4)], 3)

# Regressions with respondent gender dummy reversed (for Figure A4)
nomination.gender.results.RC.rev <- lm_robust(nomination ~ woman * I(1 - R.woman) + young, 
                                              data, subset = code != 999, 
                                              clusters = code, fixed_effects = cluster)

support.gender.results.RC.rev <- lm_robust(support ~ woman * I(1 - R.woman) + young, 
                                           data, subset = code != 999, 
                                           clusters = code, fixed_effects = cluster)

electability.gender.results.RC.rev <- lm_robust(electability ~ woman * I(1 - R.woman) + young, 
                                                data, subset = code != 999, 
                                                clusters = code, fixed_effects = cluster)

# Figure A4
png("Figure_A4.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(-0.1, 0.5), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.1, 0.5, 0.1)[-2], lty = 3, col = "lightgray")
segments(0.85, nomination.gender.results.RC.rev$conf.low[1], 
         0.85, nomination.gender.results.RC.rev$conf.high[1])
points(0.85, nomination.gender.results.RC.rev$coefficients[1], pch = 19)
text(0.85, nomination.gender.results.RC.rev$conf.high[1], "Female", pos = 3)
segments(1.15, nomination.gender.results.RC$conf.low[1], 
         1.15, nomination.gender.results.RC$conf.high[1])
points(1.15, nomination.gender.results.RC$coefficients[1], pch = 21, bg = "white")
text(1.15, nomination.gender.results.RC$conf.low[1], "Male", pos = 1)
segments(1.85, support.gender.results.RC.rev$conf.low[1], 
         1.85, support.gender.results.RC.rev$conf.high[1])
points(1.85, support.gender.results.RC.rev$coefficients[1], pch = 19)
segments(2.15, support.gender.results.RC$conf.low[1], 
         2.15, support.gender.results.RC$conf.high[1])
points(2.15, support.gender.results.RC$coefficients[1], pch = 21, bg = "white")
segments(2.85, electability.gender.results.RC.rev$conf.low[1], 
         2.85, electability.gender.results.RC.rev$conf.high[1])
points(2.85, electability.gender.results.RC.rev$coefficients[1], pch = 19)
segments(3.15, electability.gender.results.RC$conf.low[1], 
         3.15, electability.gender.results.RC$conf.high[1])
points(3.15, electability.gender.results.RC$coefficients[1], pch = 21, bg = "white")
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.gender.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.gender.summary.RC$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.gender.summary.RC$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Heterogeneity analysis by the magnitude score ####
# Note: These analyses use coarsened variables.
nomination.M.results <- lm_robust(nomination ~ woman * (magnitude + DID) + young, 
                                  data, clusters = code, 
                                  fixed_effects = cluster)
nomination.M.summary <- summary(nomination.M.results)
round(nomination.M.summary$coefficients[, c(1, 2, 4)], 3)

support.M.results <- lm_robust(support ~ woman * (magnitude + DID) + young, 
                               data, clusters = code, 
                               fixed_effects = cluster)
support.M.summary <- summary(support.M.results)
round(support.M.summary$coefficients[, c(1, 2, 4)], 3)

electability.M.results <- lm_robust(electability ~ woman * (magnitude + DID) + young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
electability.M.summary <- summary(electability.M.results)
round(electability.M.summary$coefficients[, c(1, 2, 4)], 3)

# Compute conditional ATEs by the magnitude score
median.DID <- median(data$DID, na.rm = TRUE)
nomination.M.plot <- conditional.effects.fixed(nomination.M.results, 
                                               "woman", "magnitude", 1:8, 
                                               "DID", median.DID)
support.M.plot <- conditional.effects.fixed(support.M.results, 
                                            "woman", "magnitude", 1:8, 
                                            "DID", median.DID)
electability.M.plot <- conditional.effects.fixed(electability.M.results, 
                                                 "woman", "magnitude", 1:8, 
                                                 "DID", median.DID)

# Create background histogram for Figure A7
hist.magnitude <- hist(data$magnitude, breaks = 0:8, plot = FALSE)

# Figure A11 (uses coarsened variable; not a reproduction of Figure A7)
png("Figure_A11.png", width = 6, height = 2, units = "in", pointsize = 8, res = 2000)
layout(matrix(1:3, 1, 3))
par(mar = c(3.5, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 8.5), ylim = c(-0.2, 0.4), 
     main = "Nomination willingness", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.magnitude$breaks) - 1)) {
  polygon(c(hist.magnitude$breaks[i] + 0.5, hist.magnitude$breaks[i] + 0.5, 
            hist.magnitude$breaks[i + 1] + 0.5, hist.magnitude$breaks[i + 1] + 0.5), 
          c(-0.2, hist.magnitude$density[i] / 2 - 0.2, 
            hist.magnitude$density[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:8, nomination.M.plot[, 2], 1:8, nomination.M.plot[, 3]) 
points(1:8, nomination.M.plot[, 1], pch = 19)
text(1.4, 0.37, 
     paste0("N = ", format(nomination.M.summary$nobs, big.mark = ",")))
mtext("Magnitude score", 1, 2.5, cex = 0.8)
axis(1, at = 1:8, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 8.5), ylim = c(-0.2, 0.4), 
     main = "Perceived support", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.magnitude$breaks) - 1)) {
  polygon(c(hist.magnitude$breaks[i] + 0.5, hist.magnitude$breaks[i] + 0.5, 
            hist.magnitude$breaks[i + 1] + 0.5, hist.magnitude$breaks[i + 1] + 0.5), 
          c(-0.2, hist.magnitude$density[i] / 2 - 0.2, 
            hist.magnitude$density[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:8, support.M.plot[, 2], 1:8, support.M.plot[, 3]) 
points(1:8, support.M.plot[, 1], pch = 19)
text(1.4, 0.37, 
     paste0("N = ", format(support.M.summary$nobs, big.mark = ",")))
mtext("Magnitude score", 1, 2.5, cex = 0.8)
axis(1, at = 1:8, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 8.5), ylim = c(-0.2, 0.4), 
     main = "Perceived electability", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.4, 0.1)[-3], lty = 3, col = "lightgray")
for (i in 1:(length(hist.magnitude$breaks) - 1)) {
  polygon(c(hist.magnitude$breaks[i] + 0.5, hist.magnitude$breaks[i] + 0.5, 
            hist.magnitude$breaks[i + 1] + 0.5, hist.magnitude$breaks[i + 1] + 0.5), 
          c(-0.2, hist.magnitude$density[i] / 2 - 0.2, 
            hist.magnitude$density[i] / 2 - 0.2, -0.2), 
          border = NA, col = "darkgray")
}
segments(1:8, electability.M.plot[, 2], 1:8, electability.M.plot[, 3]) 
points(1:8, electability.M.plot[, 1], pch = 19)
text(1.4, 0.37, 
     paste0("N = ", format(electability.M.summary$nobs, big.mark = ",")))
mtext("Magnitude score", 1, 2.5, cex = 0.8)
axis(1, at = 1:8, lwd = 0.5, cex = 0.8)
mtext("Effect of female candidate", 2, 2.5, cex = 0.8)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Subgroup analyses of ordinance-designated cities ####
nomination.ODC.results <- lm_robust(nomination ~ woman + young, 
                                    data, subset = designated == 1, 
                                    clusters = code, 
                                    fixed_effects = cluster)
nomination.ODC.summary <- summary(nomination.ODC.results)
round(nomination.ODC.summary$coefficients[, c(1, 2, 4)], 3)

support.ODC.results <- lm_robust(support ~ woman + young, 
                                 data, subset = designated == 1, 
                                 clusters = code, 
                                 fixed_effects = cluster)
support.ODC.summary <- summary(support.ODC.results)
round(support.ODC.summary$coefficients[, c(1, 2, 4)], 3)

electability.ODC.results <- lm_robust(electability ~ woman + young, 
                                      data, subset = designated == 1, 
                                      clusters = code, 
                                      fixed_effects = cluster)
electability.ODC.summary <- summary(electability.ODC.results)
round(electability.ODC.summary$coefficients[, c(1, 2, 4)], 3)

# Figure A8
png("Figure_A8.png", width = 3, height = 3, units = "in", pointsize = 8, res = 2000)
par(mar = c(2.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 3.5), ylim = c(-0.2, 0.6), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, col = "lightgray")
abline(h = seq(-0.2, 0.6, 0.2)[-2], lty = 3, col = "lightgray")
segments(1, nomination.ODC.results$conf.low[1], 
         1, nomination.ODC.results$conf.high[1])
points(1, nomination.ODC.results$coefficients[1], pch = 19)
segments(2, support.ODC.results$conf.low[1], 
         2, support.ODC.results$conf.high[1])
points(2, support.ODC.results$coefficients[1], pch = 19)
segments(3, electability.ODC.results$conf.low[1], 
         3, electability.ODC.results$conf.high[1])
points(3, electability.ODC.results$coefficients[1], pch = 19)
mtext(c(paste0("Nomination\nwillingness\n(N = ", 
               format(nomination.ODC.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nsupport\n(N = ", 
               format(support.ODC.summary$nobs, big.mark = ","), ")"), 
        paste0("Perceived\nelectability\n(N = ", 
               format(electability.ODC.summary$nobs, big.mark = ","), ")")), 
      1, 1.5, at = 1:3)
mtext("Average treatment effect of female candidate", 2, 2.5)
axis(2, lwd = 0.5, cex = 0.8)
dev.off()

#### Heterogeneity analysis by respondents' years in office ####
# Table A8
# Note: These analyses use coarsened variables and do not replicate Table A7.
nomination.tenure.results <- lm_robust(nomination ~ woman * (R.age + tenure) + young, 
                                       data, clusters = code, 
                                       fixed_effects = cluster)
nomination.tenure.summary <- summary(nomination.tenure.results)
round(nomination.tenure.summary$coefficients[, c(1, 2, 4)], 3)
nomination.tenure.results$nobs

support.tenure.results <- lm_robust(support ~ woman * (R.age + tenure) + young, 
                                    data, clusters = code, 
                                    fixed_effects = cluster)
support.tenure.summary <- summary(support.tenure.results)
round(support.tenure.summary$coefficients[, c(1, 2, 4)], 3)
support.tenure.results$nobs

electability.tenure.results <- lm_robust(electability ~ woman * (R.age + tenure) + young, 
                                         data, clusters = code, 
                                         fixed_effects = cluster)
electability.tenure.summary <- summary(electability.tenure.results)
round(electability.tenure.summary$coefficients[, c(1, 2, 4)], 3)
electability.tenure.results$nobs